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Abstract 

Computational time reversal imaging can be used to locate the position 
of multiple scatterers in a known background medium. The current meth- 
ods for computational time reversal imaging are based on the null subspace 
projection operator, obtained through the singular value decomposition 
of the frequency response matrix. Here, we discuss the image recovery 
problem from a small number of random and noisy measurements, and 
we show that this problem is equivalent to a randomized approximation 
of the null subspace of the frequency response matrix. 
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1 Introduction 



Computational time-reversal imaging (CTRI) has become an important research 
area in recent years, with relevant applications in radar imaging, exploration 
seismics, nondestructive material testing, medical imaging [1-9] etc. CTRI uses 
the information carried by scattered acoustic, elastic or electro-magnetic waves 
to obtain images of the investigated domain [1]. It was shown that scattered 
acoustic waves can bo time-reversed and focused onto their original source lo- 
cation through arbitrary media, using a so-called time- reversal mirror [2] . This 
important result shows how one can use CTRI to identify the location of mul- 
tiple point scatterers (targets) in a known background medium [3]. In this 
case, a back-propagated signal is computed, rather than implemented in the 
real medium, and its peaks indicate the existence of possible scattering targets. 
The current methods for CTRI are based on the null subspace projection opera- 
tor, obtained through the singular value decomposition (SVD) of the frequency 
response matrix [4-9]. Motivated by several results obtained in random low 
rank approximation theory, here we investigate the problem of image recovery 
from a small number of random and noisy measurements, and we show that this 
problem is equivalent to a randomized approximation of the null subspace of 
the frequency response matrix. 

2 Frequency response matrix 

We consider a system consisting of an array of transceivers (i.e. each antenna 
is an emitter and a receiver) located at Xn G (n = 1, N), and a collection 
of M distinct scatterers (targets) with scattering coefficients p,„. located at 
Um G {tti = 1, M) (Fig. 1). Here, _D = 1, 2, 3 is the dimensionality of the 
space. Also, we assume that the wave propagation is well approximated in the 
space-frequency domain {x,uj) by the inhomogeneous Helmholtz equation [1-9]: 

[^'^ + klrf{x)\i}{x,u) = -s{x,u:), (1) 

where i^(x, uj) is the wave amplitude produced by a localized source s(x, u), ko = 
27rw/co = 2tt/X is the wavenumber of the homogeneous background, with lu the 
frequency, Cq the homogeneous background wave speed, and A the wavelength. 
Here, rj^x) is the index of refraction: r]{x) = Co/c{x), where c{x) is the wave 
speed at location x. In the background we have 77o(a;) = 1, while 'rf'{x) = 
1 -|- a{x), measures the change in the wave speed at the scatterers location. 

The fundamental solutions, or the Green functions, for this problem satisfy 
the following equations: 

[V' + kl] Go{x, x') = -S{x - x'), (2) 

+ kWix)] G{x, x') = -5{x - x'), (3) 

for the homogeneous and inhomogeneous media, respectively. The fundamental 
solution G{x, x') for the inhomogeneous medium can be written in terms of that 
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for the homogeneous one Gq{x,x') as: 

G{x,x')=Go{x,x') + kl j a{z)Go{x,z)G{z,x')dz. (4) 

This is an impHcit integral equation for G{x,x'). Since the scatterers are as- 
sumed to be pointUke, the regions with a{z) ^ are assumed to be finite, and 
included in compact domains Q.m centered &t ym, m = , which are small 

compared to the wavelength A. Therefore we can write: 

M 

a{z,oj) = ^ Pm{<^)S{z - ym), (5) 

m=l 

and consequently we obtain: 

M 

G{x, x') ~ Go{x, x') + ^ pm{'^)Go{x, ym)G{ym, x'). (6) 

m— 1 

If the scatterers are sufficiently far apart we can neglect the multiple scattering 
among the scatterers {G{ym,x') ~ Go{ym,x')) and we obtain the Born approx- 
imation of the solution [10] : 

M 

G{x,x') Go{x,x') + ^ pmii^)Goix,ym)Go{ym,x'). (7) 

m=l 

If x corresponds to the receiver location Xi, and x' corresponds to the emitter 
location Xj , then we obtain: 

G{xi , a;j ) ~ Go {xi ,Xj)+ Hij (w) , (8) 

where 

M 

Hij{u))='^Go{Xi,ym)pm{i^)Go{ym,Xj), i,j = l,...,N, (9) 

m—l 

are the elements of the frequency response matrix H{uj) = [Hij{uj)]. The re- 
sponse matrix H(uj) is obviously a complex and symmetric N x N matrix, since 
the same Green function is used in both the transmission and the reception 
paths. 

3 Computational time-reversal imaging 

An important step in CTRI is to determine the frequency response matrix H^uj). 
This can be done by performing a series of N experiments, in which a single 
element of the array is excited with a suitable signal s{uj) and wc measure the 
frequency response between this element and all the other elements of the array 
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[1-9]. In general, given the Green function Go{x,x'), the general solution to the 
Helmholtz equation is the convolution; 

il){x,Lo) = {Go * s){x,uj) = J Go{x,x')s{x' ,w)dx' . (10) 

Thus, if the j antenna emits a signal Sj{uj) then, using the convolution theorem 
in the Fourier domain, the field produced at the location r is Go{r,Xj)sj{uj). 
If this field is incident on the m-th scattcrer, it produces the scattered field 
Go{r,ym)Pm{^)Go{r,Xj)sj. Thus, the total wave field, due to a pulse emitted 
by a single element at Xj and scattered by the M targets can be expressed as: 

M 

V'(r,a;) = ^ GQ{r,y^)pm{ui)Go{y^,Xj)sj{ui). (11) 

m— 1 

If this field is measured at the i-th antenna we obtain: 

M 

tp{Xi,U)) = ^ Go{Xi,ym)Pm{i^)Go{ym,Xj)Sj = Hij{u))Sj{u}). (12) 
m=l 

In CTRI one forms the symmetric self-adjoint matrix [1-9]: 

K{oj) = H*{uj)H{w) = H{uj)H{<jj), (13) 

where the star denotes the adjoint and the bar denotes the complex conjugate 
[H* = H, since H is symmetric) . H is the frequency- domain version of a time- 
reversed response matrix, thus if(w) corresponds to performing a scattering 
experiment, time-reversing the received signals and using them as input for a 
second scattering experiment. Therefore, time- reversal imaging relies on the 
assumption that the Green function can be always calculated. 

As long as the number of transceivers exceeds the number of scatterers, M < 
N, the matrix K{u)) is rank deficient and it has only M non-zero eigenvalues, 
with the corresponding eigenvectors Vmit^), m = 1, M. When the scatterers 
are well resolved, the eigenvectors can be back-propagated as (j(^(r, cj)?;m(a;), 
and consequently the radiated wavefields focus at target locations. Thus, each 
eigenvector can be used to locate a single scatterer. Here, g{r,co) is the Green 
function vector, which expresses the response at each array element due to a 
single pulse emitted from r: 

g{r,uj) = [ GQ{xi,r,uj) Gn{x2,r,uj) ... Go{xN,r,u!) . (14) 

The above result does not apply to the case of poorly-resolved targets. In 
this case, the eigenvectors of if (w) are linear combinations of the target Green 
function vectors g{ym,oj). Thus, back-propagating one of these eigenvectors 
generates a linear combination of wavefields, each focused on a different target 
location. The subspace-based algorithms, based on the multiple signal classifi- 
cation (MUSIC) method, can be used in this more general situation [7-9] . The 
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signal subspacc method assumes that the mimbcr M of point targets in the 
medium is lower than the number of transceivers N , and the general idea is to 
localize multiple sources by exploiting the eigenstructure and the rank deficiency 
of the response matrix Hiuo). 

The SVD of the matrix Kioj) is given by: 

N 

K{U}) ^^Xn{uj)Un{uj)v*^{u}), (15) 
n=l 

where Un{oj) and Wn(a;) are the left and right singular vectors. Since H{uj) is 
rank-deficient, all but the first M singular values vanish: Ai(a;) > ... > Am('^) > 
0, Xj{uj) = 0, j = M + 1,...,N. Therefore, the first M singular vectors span 
the essential signal subspace, while the remaining N — M columns span the 
null-subspace. The projection on the null-subspace is given by: 

N M 
Pnuu{'^)= X! ^n{'^)v*n{Lo) = I -^Vn{oj)vl{uj), (16) 
n=M+l n=l 

where / is the identity matrix. It follows immediately that Pnun{^)K{w) = 0, 

and therefore Pnuu{'^)<){f,ijj) = 0, for any u. Therefore, the target locations 
must correspond to the peaks in the MUSIC pseudo-spectrum for any oj: 

SMUsic{r,oj) = \\Pnuii{i^)9ir,uj)\\~'^ , (17) 

where g{r,uj) is the free-space Green function vector. Thus, one can form an 
image of the scatterers by plotting, at each point r, the quantity Smusic{''',(jj). 
The resulting plot will have large peaks at the locations of the scatterers. 



4 Randomized null-subspace approximation 

Let us first present several results on random rank approximation. The SVD of 
a. M X N {N < M) matrix A = [aij] can be written as [11]: 



(18) 



where R < N is the rank of A, Ur and Vr are the left and right singular vectors, 
and the singular values (in decreasing order) are: Ai > A2 > ... > Ai? > 0. If 
N > M one can always compute the SVD of the transpose matrix and then swap 
the left and right singular vectors in order to recover the SVD of the original 
matrix. We also remind that the Probenius and the spectral norms of A are 
[11]: 



M N 



m=l n=l 



11^112 = Ai. 



(19) 
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If we define 

K 



Ak = J2^kUkvl, (20) 

fc=i 

for any K < R, then, by the Eckart- Young theorem, Ak is the best ranlc K 
approximation to A with respect to the spectral norm and the Probenius norm 
[11]. Thus, for any matrix B of rank at most K, we have: 

\\A-Ak\\1<\\A-B\\1, \\A-Ak\\1<\\A-B\\1 (21) 

From basic linear algebra we have: 

Ak = A 



■ K 

.fc=l 



(22) 



Also, wc say that a matrix A has a good rank K approximation if A — Ak is 
small with respect to the spectral norm and the Frobenius norm. 

Our problem is to substitute A^ with some other rank K matrix D , which is 
much simpler than Ax, and does not require the full knowledge of A. Therefore, 
the matrix D must satisfy the general condition: 

\\A-D\\1<\\A-Ak\\1+^, (23) 

where ^ represents a tolerable level of error for the given application. Several 
important results have been recently obtained regarding this problem. 

It has been shown that one can compute a rank K approximation of A from 
a randomly chosen submatrix of A [12, 13]. For any K < R and Q < e,5 < 1 
this method uses a matrix D, containing only a random sample of K rows of 
matrix A, so that: 

\\A-D\\l<\\A-AK\\l+e\\A\\l, (24) 

holds with probability of at least 1 — 5. Recently, the above result has been 
improved: 

\\A-D\\l<{l + e)\\A-AK\\l, (25) 

by taking into account that the additive error e ||^||^ can be arbitrarily large 
compared to the true error ||^ — ^ifUfr [14]. These results show that the sparse 
matrix D recovers almost as much from A as the best rank approximation matrix 
Ak. 

In a different approach [15], it has been shown that one can substitute Ak 

with a sparse matrix D = [dij], where dij = aij with probability p, and dij = 
with probability I — p. This result asserts that it is possible to find a good low 
rank approximation to A even after randomly omitting many of its entries. In 
particular, it has been shown that the stronger the spectral features of A the 
more of its entries we can afford to omit. 

Another observation is related to the noise effect on low rank approximation. 
One can model this by adding to A a matrix F whose entries are independent 
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Gaussian random variables with mean and standard deviation a. As long 
as a is not too big, the optimal rank K approximation {A + F)k to A + F 
will approximate A nearly as well as Ak [15]. This stability of low rank ap- 
proximations with respect to Gaussian noise is well-understood, and in fact low 
rank approximations are frequently used with the explicit purpose of removing 
Gaussian noise. 

The above results can be transferred to the CTRI problem by substituting 
the matrix A with the frequency response matrix H{u)), and considering a sparse 
matrix H{u) satisfying the above conditions for the matrix D. Also, from the 
CTRI considerations, we can form the matrix: 

K{uj) = H*{ijj)H{uj), (26) 

as a substitute for K{w). Since -ff(w) is a rank M approximation of H{u)), the 
rank of K{lo) will also be M, and its SVD will be given by: 

N 

K{u) = Hn{i^)an{ivK{iv), (27) 

n=l 

where a„(ci;) and bn{Lo) are the left and right singular vectors, and all but the 
first M singular values vanish: Hi{u;) > ... > ^m{'.^) > 0, l^j{oj) = 0, j = 
M + 1, ...,N. Now, since H{ui) is a rank M approximation of H{u)), the last 
N — M right singular vectors of K{(jj) will approximate the null-subspace of 
K{uj) (and implicitly of H{w)). Thus, the approximate projection on the null- 
subspace is 

JV M 

n=M+l Ti=l 

with the corresponding MUSIC pseudo-spectrum given by: 

_ _ -2 

SMUsic{r,i^) = Pnuii{i^)g{r,u)) . (29) 

In order to illustrate and validate numerically the above results, we have con- 
sidered a two dimensional scenario, including N = 100 transceivers, separated 
hy d= A/2 and located at a;„ = [ nX/2 + a/2 - NX/ A where a = lOOA 
is the side of the imaging area. The number of targets (with the scattering 
coefficients pm = 1) is set to M = 5 and their position is randomly generated in 
the imaging area. The computational image grid is also set to L x L = 300 x 300 
pixels. The two dimensional Green function is Go{x,x') = jH'^\kQ\x — x'\), 

where Hq^\.) is the zero order Hankel function of the first kind. The noise level 
is characterized by the signal to noise ratio (SNR). SNR compares the level of 
a desired signal to the level of background noise. The higher the ratio, the less 
obtrusive the background noise is. SNR measures the power ratio between a 
signal and the background noise: 

SJVR = Psignal / Pnoise — i^Asignal / A^oise) 5 (^^) 
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where P is average power and A is root mean sqiiare (RMS) amplitude. 

Let us first consider the case when the matrix H{uj) is obtained by randomly 
selecting M < J < N rows of the matrix H{uj). In Figure 2 we give the results 
obtained for the extreme case of J = M and for different levels of noise, SNR = 
00, 80, 40, 20, 10, 5, 2, 1. One can see that even for this extreme case the results 
are actually pretty good. By increasing J the quality of the image improves 
even at high levels of noise, as shown on Figure 3, where J — IM, 2M, 3M, 4M 
and the noise is fixed at SNR = 2 on the first line of images, and respectively 
SNR = 1 on the second line of images. If J < M the algorithm doesn't work. 

In Figure 4 we give the results obtained when the elements of the ma- 
trix H{uj) are selected randomly as: Hij(u)) = Hij(u>) with probability p, and 
Hiji^bj) = with probability 1 — p (we also conserved the symmetry Hij{Lo) = 
Hji{uj) during the random selection). The figure 'matrix' is organized on lines 
and columns. The lines correspond to the probability p = 1,0.5,0.25,0.1, and 
the columns correspond to the noise level SNR = cx), 10, 5, 2. 

5 Conclusion 

We have shown that the problem of image recovery from a small number of ran- 
dom and noisy measurements is equivalent to a randomized approximation of 
the null subspace of the frequency response matrix. The obtained results show 
that one can recover the sparse time-reversal image from fewer (random) mea- 
surements than conventional methods use. Prom the analytical results and the 
numerical experiments we conclude that the minimum number of measurements 
is MN <C 7V^, where M is the rank of the full matrix H{lo). 
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Figure 1: Geometry of a time-reversal imaging experiment, containing 
transceivers and M scattering targets. 



10 



Figure 2: Numerical results for H{u!) obtained by randomly selecting 
J = M rows of the matrix H{uj), for different levels of noise: SNR = 
00, 80, 40, 20, 10, 5, 2, 1 (from top left corner to bottom right corner). 
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Figure 3: Numerical results for -ff(a;) obtained by randomly selecting J = 
M,2M,3M,4M rows of the matrix H{oj): first line SNR = 2; second line 
SNR=1. 
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Figure 4: Numerical results obtained when the elements of the matrix H{lo) are 
selected randomly as: Hij{uj) = Hij{u}) with probability p, and Hij{uj) = with 
probability 1 — p. The lines correspond to the probability p = 1,0.5,0.25,0.1, 
and the columns correspond to the noise level SNR = oo, 10, 5, 2. 
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